%% Setup
% Information Structure
Info            =   struct();
Info.ParamList  =   {'RHO_M_dom','RHO_M_ext','RHO_Y_H',...
                     'Taylor_Pi','Taylor_ir',...'Cons_ElasTN','Cons_ElasHF',...
                     'Fin_AdjCost','Cons_Elas'};
Info.ParamNum   =   length(Info.ParamList);
Info.StdVec     =   [1 1 1 -1 -1 60 -10]'*0.01/4;
Info.ValueList  =   struct('RHO_M_dom',(0.05:0.1:0.95)', ...
                           'RHO_M_ext',(0.05:0.1:0.95)', ...
                           'RHO_Y_H',(0.05:0.1:0.95)',...
                           'Taylor_Pi',(1.00:0.10:2.0)',...
                           'Taylor_ir',(0:0.1:0.9)',...'Cons_ElasTN',(1:1:10)',...'Cons_ElasHF',(1:1:10)',...
                           'Fin_AdjCost',(0.5:0.5:8)',...
                           'Cons_Elas',(1:1:10)');
Info.FixedParam =   struct('RHO_M_dom',0.68,'RHO_M_ext',0.81,'RHO_Y_H',0.95,...
                           'Taylor_Pi',1.1,'Taylor_ir',0.87,...'Cons_ElasTN',6.19,'Cons_ElasHF',6.19,...
                           'Fin_AdjCost',4.75,'Cons_Elas',6.19);
% Info.FixedParam =   struct();
Info.Num        =   0;
UnitValue       =   zeros(1,Info.ParamNum);
for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    Info.Idx.(vv)   =   Info.Num+(1:1:length(Info.ValueList.(vv)))';
    Info.Num        =   Info.Num+length(Info.ValueList.(vv));
    UnitValue(ii)   =   Info.FixedParam.(vv);
end

Info.BenchMark  =   UnitValue;
ParamMat        =   repmat(UnitValue,[Info.Num,1]);

for ii=1:Info.ParamNum
    vv              =   Info.ParamList{ii};
    ParamMat(Info.Idx.(vv),ii) ...
                    =   Info.ValueList.(vv);
end

%% Computation
TempFun         =   @(Param)SubFun_Param2IRF(PP,SS,MODEL,EquJac,Info,Param);

Result          =   cell(Info.Num,1);

parfor ii=1:Info.Num
    Result{ii}      =   TempFun(ParamMat(ii,:));
end

Result_Benchmark=   TempFun(Info.BenchMark);

save('Identification.mat','Info','Result','Result_Benchmark');

return
%% Collect the Results
ParamLabel      =   struct('RHO_M_dom','$\rho_{M}$',...
                           'RHO_Y_H','$\rho^{H}_{Y}$',...
                           'Taylor_Pi','$\lambda_{\pi}$',...
                           'Taylor_ir','$\lambda_{i}$',...
                           'Cons_Elas','$\theta$',...
                           'Fin_AdjCost','$\phi$');
StatList        =   {'C to $\varepsilon_{M}$','ir to $\varepsilon_{M}$',...
                     'UIP Dev to $\varepsilon_{M}$','Rel Price to $\varepsilon_{M}$',...
                     'Exchange Rate/C to $ \varepsilon^{H}_{Y} $',...
                     'Export/C to $ \varepsilon^{H}_{Y} $'};
TempFun         =   @(IRF)[SubFun_MaxAbsVal(IRF.Eps_Eta_M_dom.C)*100;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_M_dom.ir)*400;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_M_dom.DevUIP)*400;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_M_dom.RelPrice)*100;
                           SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.ExchangeRate)/SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.C);...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.Export)/SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.C)];
StatMat         =   NaN(6,Info.Num);
for ii=1:Info.Num
    if ~isempty(Result{ii})
        StatMat(:,ii)   =   TempFun(Result{ii});
    end
    
end 
BenchMarkStat   =   TempFun(Result_Benchmark);
%% Figures
Folder          =   'TableGraphs/Identification/';
mkdir(Folder);
ParamList       =   fieldnames(ParamLabel);
N_col           =   length(ParamList);
N_row           =   size(StatMat,1);
fig_size        =   [N_col,N_row].*[1,1]*1/4; 
fig_gap         =   [0.03,0.03];
fig_VMargin     =   [0.05,0.05];
fig_HMargin     =   [0.05,0.05]; 

fig_SubPlot     =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                         fig_gap,fig_VMargin,fig_HMargin);
Fig             =   FigureSetup(['Identification'],fig_size);
for ii=1:N_col
    vv          =   ParamList{ii};
    for jj=1:N_row
        ax          =   fig_SubPlot((jj-1)*N_col+ii);
        XAxis       =   Info.ValueList.(vv)';
        YAxis       =   StatMat(jj,Info.Idx.(vv));
        LineStyle   =   {'-';'b';2};
        plot(XAxis,YAxis,'-b','LineWidth',2);
        yline(BenchMarkStat(jj),':k','LineWidth',1);
        xline(Info.FixedParam.(vv),':k','LineWidth',1);
        ylabel(StatList{jj},'Interpreter','Latex','FontSize',8); 
        title(ParamLabel.(vv),'Interpreter','Latex','FontSize',8);
        xlim([min(XAxis),max(XAxis)]);
        ylim([min(min(StatMat(jj,:)),BenchMarkStat(jj)),max(max(StatMat(jj,:)),BenchMarkStat(jj))]);
        ax.FontSize =   8;
        ax.XGrid    =   'on'; ax.YGrid  =   'on';
%         ax.XAxis.Exponent   =   0;
%       ax.YAxis.Exponent   =   0;
%         ax.YAxis.TickLabelFormat    =   '%g2';
    end
end
print('-depsc','-r1000',Fig,[Folder,'Identification']);


%% Wage Response Difference
% Data
ParamLabel      =   struct('RHO_M_dom','$\rho_{M}$',...
                           'RHO_M_ext','$\rho^{*}_{M}$',...
                           'RHO_Y_H','$\rho^{H}_{Y}$',...
                           'Taylor_Pi','$\lambda_{\pi}$',...
                           'Taylor_ir','$\lambda_{i}$',...
                           'Cons_ElasTN','$\theta_{TN}$',...
                           'Cons_ElasHF','$\theta_{HF}$',...
                           'Fin_AdjCost','$\phi$');
StatList        =   {'$ w^{H} $ to $ \varepsilon^{H}_{Y} $',...
                     '$ w^{N} $ to $ \varepsilon^{H}_{Y} $',...
                     '$ w^{H} -w^{N} $ to $ \varepsilon^{H}_{Y} $'};
TempFun         =   @(IRF)[SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.w_H)*100;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.w_N)*100;...
                           SubFun_MaxAbsVal(IRF.Eps_Eta_Y_H.w_H-IRF.Eps_Eta_Y_H.w_N)*100 ...
                           ];
StatMat         =   NaN(3,Info.Num);
for ii=1:Info.Num
    if ~isempty(Result{ii})
        StatMat(:,ii)   =   TempFun(Result{ii});
    end
    
end 
BenchMarkStat   =   TempFun(Result_Benchmark);

% Figure
Folder          =   'TableGraphs/Identification/';
mkdir(Folder);
ParamList       =   fieldnames(ParamLabel);
N_col           =   length(ParamList);
N_row           =   size(StatMat,1);
fig_size        =   [N_col,N_row].*[1,1]*1/4; 
fig_gap         =   [0.05,0.04];
fig_VMargin     =   [0.05,0.05];
fig_HMargin     =   [0.05,0.05]; 

fig_SubPlot     =   @(ii_vv)subtightplot(N_row,N_col,ii_vv,...
                                         fig_gap,fig_VMargin,fig_HMargin);
Fig             =   FigureSetup(['Identification'],fig_size);
for ii=1:N_col
    vv          =   ParamList{ii};
    for jj=1:N_row
        ax          =   fig_SubPlot((jj-1)*N_col+ii);
        XAxis       =   Info.ValueList.(vv)';
        YAxis       =   StatMat(jj,Info.Idx.(vv));
        LineStyle   =   {'-';'b';2};
        plot(XAxis,YAxis,'-b','LineWidth',2);
        yline(BenchMarkStat(jj),':k','LineWidth',1);
        xline(Info.FixedParam.(vv),':k','LineWidth',1);
        ylabel(StatList{jj},'Interpreter','Latex','FontSize',8); 
        title(ParamLabel.(vv),'Interpreter','Latex','FontSize',8);
        xlim([min(XAxis),max(XAxis)]);
        ylim([min(min(StatMat(jj,:)),BenchMarkStat(jj)),max(max(StatMat(jj,:)),BenchMarkStat(jj))]);
        ax.FontSize =   8;
        ax.XGrid    =   'on'; ax.YGrid  =   'on';
%         ax.XAxis.Exponent   =   0;
%       ax.YAxis.Exponent   =   0;
%         ax.YAxis.TickLabelFormat    =   '%g2';
    end
end
print('-depsc','-r1000',Fig,[Folder,'Identification_WageDiff']);
